Antje Hofmann: Learning Syntactic Relations with the Mole Task

RePsychLing in SMLP2022

1 Background

Children’s (age: 4-8 years)reaction times in a task teaching them syntactic relations.

1.1 Overview

• Original analysis is by Antje Hofmann.
• MixedModels.jl version
• Addition of new chunks illustrate
  • selection of parsimonious LMM using random-effects PCA
  • plotting conditional means
  • illustration of borrowing strength

2 Readme

2.1 Design

2.2 Variables

• Subj: Participant ID (renamed from ID; random factor)
• Item: Word ID (random factor)
• age: 4 - 8 years
• Block (within-Subj/within-Item):
  • 1st Learning
  • 2nd Learning
  • Disruption
  • Recovery
• Target(renamend fom targetness)
  • non-syllable target
  • syllable target
• rt: response time

3 Setup

3.1 Environment

3.2 Packages

First attach the MixedModels.jl package and other packages for plotting. The CairoMakie.jl package allows the Makie graphics system [@Danisch2021] to generate high quality static images. Activate that package with the SVG (Scalable Vector Graphics) backend.

Code

using AlgebraOfGraphics
using Arrow
using CairoMakie       # graphics back-end
using CategoricalArrays
using Chain
using DataFrames
using DataFrameMacros  # simplified dplyr-like data wrangling 
using MixedModels
using MixedModelsMakie # diagnostic plots
using ProgressMeter
using Random           # random number generators
using RCall            # call R from Julia
using StatsModels

using AlgebraOfGraphics: boxplot
using AlgebraOfGraphics: density

using MixedModelsMakie: qqnorm
using MixedModelsMakie: ridgeplot
using MixedModelsMakie: scatter
using MixedModelsMakie: caterpillar

ProgressMeter.ijulia_behavior(:clear);
CairoMakie.activate!(; type="svg");

• The data are available as an arrow file.
• Most of preprocessing was done with R in RStudio (see Hofmann_Maulwurf.Rmd).
• Order of factor levels should be checked.

Code

dat = DataFrame(Arrow.Table("./data/Hofmann_Maulwurf_rt.arrow"))
transform!(dat, 
    :Target => categorical => :Target,
    :Block  => categorical => :Block,
    :age => (x -> x .- 6) => :a1, # center age (linear) at 6 years
    :rt => (x -> log.(x)) => :lrt)
describe(dat)

8 rows × 7 columns

	variable	mean	min	median	max	nmissing	eltype
	Symbol	Union…	Any	Union…	Any	Int64	Type
1	Subj		p01		p85	0	Union{Missing, String}
2	Item		01PA.ogg		32SG.ogg	0	Union{Missing, String}
3	age	6.21081	3.7	6.31	8.04	0	Union{Missing, Float64}
4	Block		1st Learning		Recovery	0	Union{Missing, CategoricalValue{String, UInt32}}
5	Target		Non-target syllable		Target syllable	0	Union{Missing, CategoricalValue{String, UInt32}}
6	rt	2335.39	228.0	2100.0	8389.0	0	Union{Missing, Float64}
7	a1	0.210813	-2.3	0.31	2.04	0	Float64
8	lrt	7.67367	5.42935	7.64969	9.03468	0	Float64

• Centering age at six years yields an interpretable GM
  
• Factor levels can also be set when contrasts are defined (see below).
• BoxCox check showed that reaction time rt [ms] should be transformed to speed [1/s] = [Hz]
• Indicator variables for Target and Block generated in R.

4 LMM analysis

4.1 Contrasts

Code

contrasts = merge(
      Dict(:Target => EffectsCoding()),
      Dict(:Block => SeqDiffCoding()),
      Dict(nm => Grouping() for nm in (:Subj, :Item))
   );

4.2 Varying only intercepts LMM m_voi

Code

f_voi1 = @formula(lrt ~  1 + Block * Target * a1 + (1 | Subj) + (1 | Item));
m_voi1 = fit(MixedModel, f_voi1, dat; contrasts)

Minimizing 53    Time: 0:00:00 (15.47 ms/it)

	Est.	SE	z	p	σ_Item	σ_Subj
(Intercept)	7.6895	0.0371	207.05	<1e-99	0.0239	0.2637
Block: 2nd Learning	-0.0920	0.0104	-8.87	<1e-18
Block: Disruption	0.0516	0.0126	4.09	<1e-04
Block: Recovery	-0.0477	0.0125	-3.81	0.0001
Target: Target syllable	-0.0108	0.0040	-2.69	0.0072
a1	-0.0896	0.0324	-2.77	0.0056
Block: 2nd Learning & Target: Target syllable	-0.0057	0.0104	-0.55	0.5854
Block: Disruption & Target: Target syllable	-0.0262	0.0124	-2.11	0.0347
Block: Recovery & Target: Target syllable	0.0020	0.0123	0.16	0.8741
Block: 2nd Learning & a1	0.0347	0.0091	3.83	0.0001
Block: Disruption & a1	0.0214	0.0108	1.97	0.0486
Block: Recovery & a1	-0.0159	0.0108	-1.47	0.1409
Target: Target syllable & a1	-0.0052	0.0035	-1.46	0.1430
Block: 2nd Learning & Target: Target syllable & a1	-0.0113	0.0090	-1.25	0.2111
Block: Disruption & Target: Target syllable & a1	0.0113	0.0108	1.05	0.2933
Block: Recovery & Target: Target syllable & a1	0.0086	0.0107	0.80	0.4248
Residual	0.3004

4.3 Extract indicator variables

Code

X = modelmatrix(m_voi1)
dat.trng = X[:,2];
dat.drpt = X[:,3];
dat.rcvr = X[:,4];
dat.trgt = X[:,5];
describe(dat)

12 rows × 7 columns

	variable	mean	min	median	max	nmissing	eltype
	Symbol	Union…	Any	Union…	Any	Int64	Type
1	Subj		p01		p85	0	Union{Missing, String}
2	Item		01PA.ogg		32SG.ogg	0	Union{Missing, String}
3	age	6.21081	3.7	6.31	8.04	0	Union{Missing, Float64}
4	Block		1st Learning		Recovery	0	Union{Missing, CategoricalValue{String, UInt32}}
5	Target		Non-target syllable		Target syllable	0	Union{Missing, CategoricalValue{String, UInt32}}
6	rt	2335.39	228.0	2100.0	8389.0	0	Union{Missing, Float64}
7	a1	0.210813	-2.3	0.31	2.04	0	Float64
8	lrt	7.67367	5.42935	7.64969	9.03468	0	Float64
9	trng	-0.014562	-0.75	0.25	0.25	0	Float64
10	drpt	-0.0528019	-0.5	-0.5	0.5	0	Float64
11	rcvr	0.0524467	-0.25	-0.25	0.75	0	Float64
12	trgt	0.00331492	-1.0	1.0	1.0	0	Float64

Switch to indicator variables and refit the model.

Code

f_voi2 = @formula(lrt ~  1 + (trng+drpt+rcvr) * trgt * a1 + (1 | Subj) + (1 | Item));
m_voi2 = fit(MixedModel, f_voi2, dat; contrasts)

	Est.	SE	z	p	σ_Item	σ_Subj
(Intercept)	7.6895	0.0371	207.05	<1e-99	0.0239	0.2637
trng	-0.0920	0.0104	-8.87	<1e-18
drpt	0.0516	0.0126	4.09	<1e-04
rcvr	-0.0477	0.0125	-3.81	0.0001
trgt	-0.0108	0.0040	-2.69	0.0072
a1	-0.0896	0.0324	-2.77	0.0056
trng & trgt	-0.0057	0.0104	-0.55	0.5854
drpt & trgt	-0.0262	0.0124	-2.11	0.0347
rcvr & trgt	0.0020	0.0123	0.16	0.8741
trng & a1	0.0347	0.0091	3.83	0.0001
drpt & a1	0.0214	0.0108	1.97	0.0486
rcvr & a1	-0.0159	0.0108	-1.47	0.1409
trgt & a1	-0.0052	0.0035	-1.46	0.1430
trng & trgt & a1	-0.0113	0.0090	-1.25	0.2111
drpt & trgt & a1	0.0113	0.0108	1.05	0.2933
rcvr & trgt & a1	0.0086	0.0107	0.80	0.4248
Residual	0.3004

They are equivalent.

4.4 A zero-correlation parameter LMM m_zcp

Code

f_zcp1 = @formula(lrt ~  1 +  Block * Target * a1 + zerocorr(1 + Block * Target | Subj) + (1 + a1 | Item));
m_zcp1 = fit(MixedModel, f_zcp1, dat; contrasts);

show(issingular(m_zcp1))
VarCorr(m_zcp1)

Minimizing 365   Time: 0:00:00 ( 2.40 ms/it)
  objective:  2393.633425634436

true

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0742531	0.2724942
	Block: 2nd Learning	0.0362822	0.1904790	.
	Block: Disruption	0.0156401	0.1250605	.	.
	Block: Recovery	0.0133038	0.1153419	.	.	.
	Target: Target syllable	0.0002450	0.0156539	.	.	.	.
	Block: 2nd Learning & Target: Target syllable	0.0000000	0.0000000	.	.	.	.	.
	Block: Disruption & Target: Target syllable	0.0032559	0.0570605	.	.	.	.	.	.
	Block: Recovery & Target: Target syllable	0.0006594	0.0256785	.	.	.	.	.	.	.
Item	(Intercept)	0.0009502	0.0308250
	a1	0.0006752	0.0259853	-0.83
Residual		0.0767542	0.2770456

Again, check the equivalence.

Code

f_zcp2 = @formula(lrt ~  1 + (trng+drpt+rcvr) * trgt * a1  + 
                zerocorr(1 + (trng+drpt+rcvr) * trgt | Subj) + (1 + a1 | Item));
m_zcp2 = fit(MixedModel, f_zcp2, dat; contrasts);

show(issingular(m_zcp2))
VarCorr(m_zcp2)

Minimizing 365   Time: 0:00:00 ( 1.51 ms/it)
  objective:  2393.633425634436

true

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0742531	0.2724942
	trng	0.0362822	0.1904790	.
	drpt	0.0156401	0.1250605	.	.
	rcvr	0.0133038	0.1153419	.	.	.
	trgt	0.0002450	0.0156539	.	.	.	.
	trng & trgt	0.0000000	0.0000000	.	.	.	.	.
	drpt & trgt	0.0032559	0.0570605	.	.	.	.	.	.
	rcvr & trgt	0.0006594	0.0256785	.	.	.	.	.	.	.
Item	(Intercept)	0.0009502	0.0308250
	a1	0.0006752	0.0259853	-0.83
Residual		0.0767542	0.2770456

4.5 A complex parameter LMM m_cpx

Code

m_cpx = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  + 
                          (1 + trgt * (trng+drpt+rcvr) | Subj) + (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_cpx))    # not ok
show(m_cpx.PCA[:Subj])     # not ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_cpx)) 
VarCorr(m_cpx)

Minimizing 1328      Time: 0:00:02 ( 1.63 ms/it)

true

Principal components based on correlation matrix

 (Intercept)   1.0     .      .      .      .      .      .      .
 trgt         -1.0    1.0     .      .      .      .      .      .
 trng          0.35  -0.35   1.0     .      .      .      .      .
 drpt          0.06  -0.06   0.35   1.0     .      .      .      .
 rcvr         -0.42   0.42  -0.23  -0.48   1.0     .      .      .
 trgt & trng   0.26  -0.26  -0.12   0.09  -0.16   1.0     .      .
 trgt & drpt  -0.4    0.4   -0.07   0.14  -0.01  -0.76   1.0     .
 trgt & rcvr   0.61  -0.61   0.16   0.04  -0.3    0.48  -0.85   1.0

Normalized cumulative variances:

[0.441, 0.667, 0.8072, 0.906, 0.9581, 1.0, 1.0, 1.0]

Component loadings
 

               PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8
 (Intercept)  -0.46   0.14   0.36  -0.2   -0.28   0.12   0.05  -0.71
 trgt          0.46  -0.14  -0.36   0.2    0.28  -0.12  -0.05  -0.71
 trng         -0.19   0.43   0.17   0.76  -0.0   -0.4   -0.13  -0.0
 drpt         -0.1    0.49  -0.59   0.15  -0.19   0.57   0.13   0.0
 rcvr          0.26  -0.42   0.31   0.47  -0.39   0.51  -0.15  -0.0
 trgt & trng  -0.31  -0.36  -0.49  -0.01  -0.52  -0.34  -0.38  -0.0
 trgt & drpt   0.39   0.44   0.14  -0.29  -0.16   0.02  -0.72  -0.0
 trgt & rcvr  -0.46  -0.18  -0.04   0.08   0.6    0.33  -0.53  -0.0

Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
──────────────────────────────────────────────────
     model-dof  -2 logLik       χ²  χ²-dof  P(>χ²)
──────────────────────────────────────────────────
[1]         28  2393.6334                         
[2]         56  2353.2482  40.3853      28  0.0611
──────────────────────────────────────────────────

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0757935	0.2753062
	trgt	0.0001016	0.0100809	-1.00
	trng	0.0334971	0.1830221	+0.35	-0.35
	drpt	0.0198200	0.1407834	+0.06	-0.06	+0.35
	rcvr	0.0191718	0.1384624	-0.42	+0.42	-0.23	-0.48
	trgt & trng	0.0020619	0.0454079	+0.26	-0.26	-0.12	+0.09	-0.16
	trgt & drpt	0.0145508	0.1206268	-0.40	+0.40	-0.07	+0.14	-0.01	-0.76
	trgt & rcvr	0.0085501	0.0924670	+0.61	-0.61	+0.16	+0.04	-0.30	+0.48	-0.85
Item	(Intercept)	0.0000000	0.0000000
	a1	0.0004752	0.0217982	+NaN
Residual		0.0768486	0.2772158

The deviance improves, but we end up with an overparameterized LMM. Remove the two small VCs of the three Block x Target interactions (i.e., the third contrast of the Block factor).

4.6 A parsimonious parameter LMM m_prm

We remove one of the VC for trgt * trng contrast interaction, that is one of the three interaction terms.

Code

m_prm1 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  + 
                          (1 + trgt *(drpt + rcvr) + trng  | Subj) + 
                          (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm1))  # not ok
show(m_prm1.PCA[:Subj])   # not ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm1, m_cpx)) 
VarCorr(m_prm1)

Minimizing 1329      Time: 0:00:02 ( 1.84 ms/it)

true
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .      .
 trgt         -1.0    1.0     .      .      .      .      .
 drpt          0.04  -0.04   1.0     .      .      .      .
 rcvr         -0.41   0.41  -0.47   1.0     .      .      .
 trgt & drpt  -0.43   0.43   0.22  -0.04   1.0     .      .
 trgt & rcvr   0.63  -0.63   0.02  -0.31  -0.88   1.0     .
 trng          0.34  -0.34   0.33  -0.23  -0.12   0.18   1.0

Normalized cumulative variances:
[0.4716, 0.7125, 0.8294, 0.9405, 0.9938, 1.0, 1.0]

Component loadings
                PC1    PC2    PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.5    0.01  -0.27   0.34   0.24   0.07  -0.71
 trgt          0.5   -0.01   0.27  -0.34  -0.24  -0.07  -0.71
 drpt         -0.08   0.63   0.22  -0.41   0.61   0.09  -0.0
 rcvr          0.28  -0.48   0.39   0.32   0.64  -0.18   0.0
 trgt & drpt   0.36   0.45  -0.35   0.35   0.04  -0.65   0.0
 trgt & rcvr

  -0.47  -0.23   0.12  -0.42  -0.04  -0.73   0.0
 trng         -0.24   0.34   0.72   0.44  -0.33  -0.07   0.0

Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + drpt + rcvr + trgt & drpt + trgt & rcvr + trng | Subj) + (1 + a1 | Item)
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
───────────────────────────────────────────────────
     model-dof  -2 logLik        χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         28  2393.6334                          
[2]         48  2338.3644   55.2691      20  <1e-04
[3]         56  2353.2482  -14.8838       8     NaN
───────────────────────────────────────────────────

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0759040	0.2755068
	trgt	0.0000989	0.0099437	-1.00
	drpt	0.0193595	0.1391386	+0.04	-0.04
	rcvr	0.0188340	0.1372370	-0.41	+0.41	-0.47
	trgt & drpt	0.0106015	0.1029636	-0.43	+0.43	+0.22	-0.04
	trgt & rcvr	0.0081492	0.0902726	+0.63	-0.63	+0.02	-0.31	-0.88
	trng	0.0349529	0.1869571	+0.34	-0.34	+0.33	-0.23	-0.12	+0.18
Item	(Intercept)	0.0008589	0.0293075
	a1	0.0006867	0.0262043	-0.91
Residual		0.0763357	0.2762891

We don’t lose goodness of fit, but are still overparameterized. We remove CPs for Target.

Code

m_prm2 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  + 
                          (1 + drpt + rcvr + trng + drpt&trgt + rcvr&trgt  | Subj) + 
                  zerocorr(0 + trgt  | Subj) +
                          (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm2))  #  ok
show(m_prm2.PCA[:Subj])   #  ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm2,  m_cpx)) 
VarCorr(m_prm2)

Minimizing 1001      Time: 0:00:02 ( 2.01 ms/it)

false
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .     .
 drpt          0.04   1.0     .      .      .      .     .
 rcvr         -0.41  -0.47   1.0     .      .      .     .
 trng          0.33   0.34  -0.22   1.0     .      .     .
 drpt & trgt  -0.36   0.26  -0.1   -0.07   1.0     .     .
 rcvr & trgt   0.59  -0.01  -0.27   0.12  -0.87   1.0    .
 trgt          0.0    0.0    0.0    0.0    0.0    0.0   1.0

Normalized cumulative variances:
[0.3441, 0.5974, 0.7402, 0.8568, 0.9478, 0.9939, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.52   0.07  0.0   -0.04   0.61  -0.58   0.15
 drpt         -0.09   0.61  0.0   -0.05  -0.62  -0.47   0.11
 rcvr          0.3   -0.49  0.0    0.52  -0.16  -0.58  -0.17
 trng         -0.26   0.36  0.0    0.83   0.09   0.3   -0.08
 drpt & trgt   0.47   0.44  0.0   -0.1    0.37  -0.14  -0.65
 rcvr & trgt  -0.59  -0.22  0.0   -0.14  -0.27   0.03  -0.72
 trgt          0.0    0.0   1.0    0.0    0.0    0.0    0.0

Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + (1 + a1 | Item)
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
───────────────────────────────────────────────────
     model-dof  -2 logLik        χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         28  2393.6334                          
[2]         42  2341.3448   52.2886      14  <1e-05
[3]         56  2353.2482  -11.9034      14     NaN
───────────────────────────────────────────────────

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0757927	0.2753047
	drpt	0.0193698	0.1391753	+0.04
	rcvr	0.0188467	0.1372834	-0.41	-0.47
	trng	0.0351278	0.1874242	+0.33	+0.34	-0.22
	drpt & trgt	0.0099421	0.0997100	-0.36	+0.26	-0.10	-0.07
	rcvr & trgt	0.0066351	0.0814561	+0.59	-0.01	-0.27	+0.12	-0.87
	trgt	0.0002651	0.0162812	.	.	.	.	.	.
Item	(Intercept)	0.0008262	0.0287433
	a1	0.0006740	0.0259607	-0.93
Residual		0.0762062	0.2760546

Check the Item-related CP. It is very large, might be spurious.

Code

m_prm3 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  + 
                          (1 + drpt + rcvr + trng + drpt&trgt + rcvr&trgt  | Subj) + 
                  zerocorr(0 + trgt  | Subj) +
                  zerocorr(1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm3))  # ok
show(m_prm3.PCA[:Subj])   # ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm3, m_prm2, m_cpx)) 
VarCorr(m_prm3)

Minimizing 1129      Time: 0:00:02 ( 2.00 ms/it)

false
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .     .
 drpt          0.04   1.0     .      .      .      .     .
 rcvr         -0.41  -0.47   1.0     .      .      .     .
 trng          0.33   0.34  -0.21   1.0     .      .     .
 drpt & trgt  -0.37   0.25  -0.08  -0.07   1.0     .     .
 rcvr & trgt   0.61   0.04  -0.31   0.13  -0.86   1.0    .
 trgt          0.0    0.0    0.0    0.0    0.0    0.0   1.0

Normalized cumulative variances:
[0.35, 0.5979, 0.7408, 0.8585, 0.9488, 0.9943, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.52   0.03  0.0   -0.01   0.61  -0.57   0.16
 drpt         -0.12   0.62  0.0   -0.06  -0.61  -0.47   0.12
 rcvr          0.33  -0.47  0.0    0.52  -0.19  -0.58  -0.18
 trng         -0.27   0.36  0.0    0.83   0.08   0.31  -0.08
 drpt & trgt   0.45   0.47  0.0   -0.1    0.38  -0.15  -0.63
 rcvr & trgt  -0.58  -0.22  0.0   -0.15  -0.26   0.03  -0.72
 trgt          0.0    0.0   1.0    0.0    0.0    0.0    0.0

Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + MixedModels.ZeroCorr((1 + a1 | Item))
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + (1 + a1 | Item)
4: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
───────────────────────────────────────────────────
     model-dof  -2 logLik        χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         28  2393.6334                          
[2]         41  2356.2006   37.4328      13  0.0004
[3]         42  2341.3448   14.8558       1  0.0001
[4]         56  2353.2482  -11.9034      14     NaN
───────────────────────────────────────────────────

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0760499	0.2757714
	drpt	0.0193875	0.1392389	+0.04
	rcvr	0.0189413	0.1376274	-0.41	-0.47
	trng	0.0350733	0.1872787	+0.33	+0.34	-0.21
	drpt & trgt	0.0091677	0.0957481	-0.37	+0.25	-0.08	-0.07
	rcvr & trgt	0.0059753	0.0773002	+0.61	+0.04	-0.31	+0.13	-0.86
	trgt	0.0002813	0.0167715	.	.	.	.	.	.
Item	(Intercept)	0.0006681	0.0258480
	a1	0.0005356	0.0231424	.
Residual		0.0763232	0.2762665

Perhaps not. We stay with LMM m_prm2

4.7 Compare models with goodness-of-fit statistics.

Code

let mods = [m_zcp2, m_prm2, m_prm1, m_cpx];
 DataFrame(;
    model=[:m_zcp2, :m_prm2, :m_prm1, :m_cpx],
    pars=dof.(mods),
    geomdof=round.(Int, (sum ∘ leverage).(mods)),
    AIC=round.(Int, aic.(mods)),
    AICc=round.(Int, aicc.(mods)),
    BIC=round.(Int, bic.(mods)),
  )
end

4 rows × 6 columns

	model	pars	geomdof	AIC	AICc	BIC
	Symbol	Int64	Int64	Int64	Int64	Int64
1	m_zcp2	28	290	2450	2450	2639
2	m_prm2	42	287	2425	2426	2709
3	m_prm1	48	276	2434	2435	2759
4	m_cpx	56	259	2465	2466	2843

LMM m_prm2 is a defensible solution according to \(\Delta\) AIC, \(\Delta\) BIC suggests we should not bother with CPs.

Code

coeftable(m_prm2)

	Coef.	Std. Error	z	Pr(>
(Intercept)	7.68255	0.0387548	198.23	<1e-99
trgt	-0.00948222	0.00442195	-2.14	0.0320
trng	-0.0841465	0.0279289	-3.01	0.0026
drpt	0.0428818	0.0228927	1.87	0.0610
rcvr	-0.0370297	0.0226366	-1.64	0.1019
a1	-0.0922379	0.0338516	-2.72	0.0064
trgt & trng	-0.00176262	0.00963315	-0.18	0.8548
trgt & drpt	-0.019929	0.01816	-1.10	0.2725
trgt & rcvr	-0.00411596	0.0162158	-0.25	0.7996
trgt & a1	-0.00874271	0.0038894	-2.25	0.0246
trng & a1	0.0414053	0.0244065	1.70	0.0898
drpt & a1	0.0287421	0.0199004	1.44	0.1487
rcvr & a1	-0.0208157	0.0197058	-1.06	0.2908
trgt & trng & a1	-0.00234843	0.00840933	-0.28	0.7800
trgt & drpt & a1	0.0171553	0.0158677	1.08	0.2796
trgt & rcvr & a1	0.00435249	0.0141867	0.31	0.7590

4.8 Pruning fixed-effects

There is an option of pruning some higher-order interactions.

4.9 Diagnostic plots

Various visualizations are used to check whether or not data are defensibly modeled with an LMM. They may lead to removal of outliers, transformations of the dependent variable, and deliver valuable heuristic information to be followed up with exploratory post-hoc analyses or ideally replication of new insights gained this way. In practice, it appears that only severe violations will stop people from reporting a model.

4.9.1 Residuals over fitted

Code

scatter(fitted(m_prm2), residuals(m_prm2))

Looks like we missed some fast response times.

4.9.2 Q-Q plot

Code

qqnorm(m_prm2; qqline=:none)

Hm?

4.9.3 Residual distributions: observed vs. theoretical

Curves for residulas based on observed and theoretical values should correspond.

Code

let
  n = nrow(dat)
  dat_rz = (;
    value=vcat(residuals(m_prm2) ./ std(residuals(m_prm2)), randn(n)),
    curve=repeat(["residual", "normal"]; inner=n),
  )
  draw(
    data(dat_rz) *
    mapping(:value; color=:curve) *
    density(; bandwidth=0.1);
  )
end

Figure 1: Kernel density plot of the standardized residuals for model m1 versus a standard normal

They are a bit too narrow.

4.10 Conditional means of random effects

4.10.1 Subject-related conditional means of random effects

Code

cm_Subj = first(ranefinfo(m_prm2))
caterpillar!(Figure(; resolution=(800, 1200)), cm_Subj; orderby=1)

Figure 2: Prediction intervals on subject random effects for model m_prm2

4.10.2 Borrowing-strength plots

Shrinkage refers to the adjustment of subject-level or item-level predictions by taking population estimates into account. The further a subject’s/item’s estimate is from the fixed effect or the more variable or less reliable the subject’s/item’s estimate, the more the prediction will be shrunk towards the population estimate. Alternative terms for shrinkage are “borrowing strength” (Tukey) and regularization. My favorite is actually Tukey’s because indeed we borrow strength from the population estimates to make predictions for individual subjects’ effects. The goal of this section to illustrate the results of borrowing strength.

Subject-related conditional means of random effects revealed information about individual differences beyond fixed effects. Would these results also be visible in unconditional means, that is when we compute GM and experimental effects within subjects (i.e., as fixed effects) without borrowing strength from the population estimates?

In the following plots, effect estimates based on alone on each subject’s data (i.e., no pooling of data, no borrowing of strength) are plotted in pink and the subjects’ conditional means shown in the caterpillar plots are plotted in blue. The arrows indicate how much a subject’s prediction is changed by borrowing strength from knowledge of the population estimates.

Code

shrinkageplot!(Figure(; resolution=(1000, 1200)), m_prm2)

Figure 3: Shrinkage plots of the subject random effects in model m_prm2

There is an outlier subject

Code

cm = raneftables(m_prm2);   
sort(DataFrame(cm.Subj), ["(Intercept)", :trng, "drpt & trgt"])
sort(DataFrame(cm.Subj), :trng)

52 rows × 8 columns (omitted printing of 1 columns)

	Subj	(Intercept)	drpt	rcvr	trng	drpt & trgt	rcvr & trgt
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	p10	-0.665206	-0.0800224	0.176762	-0.502942	0.325148	-0.281849
2	p14	-0.517114	0.14949	-0.104271	-0.380446	0.0865683	-0.0397183
3	p76	0.0792766	-0.100591	0.104864	-0.309562	-0.0285042	0.0205178
4	p02	-0.621993	-0.196999	0.260501	-0.251672	0.0704512	-0.126543
5	p36	0.421384	-0.0365387	0.0872904	-0.231847	-0.169042	0.145319
6	p39	0.291843	-0.130884	0.143621	-0.200882	-0.0192362	0.00243074
7	p74	-0.32949	-0.0607274	0.091318	-0.179369	0.0671451	-0.0858905
8	p67	0.0383213	0.00115633	-0.0194392	-0.155016	0.0331431	-0.0117934
9	p06	-0.237699	0.139392	0.026263	-0.117284	-0.0667953	0.0402144
10	p54	0.313189	-0.162019	0.0466067	-0.116629	-0.0367766	0.0411771
11	p23	-0.0288494	-0.0711022	0.0526879	-0.103847	-0.127008	0.0803724
12	p42	-0.0973896	0.288231	-0.166534	-0.10372	0.121435	-0.0338403
13	p59	0.421663	0.0211065	-0.255664	-0.10291	-0.0123478	0.0649439
14	p45	-0.363685	-0.295754	0.169737	-0.097615	0.0435265	-0.0933103
15	p30	-0.329959	-0.149979	-0.0897719	-0.0879097	-0.0134295	-0.00936962
16	p83	-0.230431	-0.0531665	-0.160467	-0.0702652	0.0132008	-0.00841336
17	p84	-0.389883	-0.202825	0.444848	-0.0491713	-0.199372	0.0205808
18	p72	0.176311	-0.0685062	-0.0523226	-0.0387344	-0.119497	0.118113
19	p12	-0.148148	-0.0448302	0.0523053	-0.0284382	-0.0671719	0.0361489
20	p04	0.138247	-0.0448374	0.0391417	-0.0278217	-0.00323117	0.00442783
21	p77	0.134277	-0.0965476	0.0522975	-0.0217883	-0.0801099	0.0500009
22	p79	-0.0769189	-0.0192919	0.0300986	-0.0144482	-0.0321694	0.010373
23	p69	0.104371	-0.0392821	-0.0829739	-0.0144277	0.0226718	0.0351996
24	p55	-0.198404	-0.0636382	-0.0785629	-0.00935968	-0.0624377	0.024083
25	p27	0.0253214	-0.00349694	-0.00408308	-0.00611823	-0.0976617	0.0736573
26	p28	0.150439	-0.00601184	-0.151526	0.00703225	0.0683167	-0.0248102
27	p25	0.267625	0.137029	-0.0474469	0.00970517	-0.0356727	0.0451016
28	p64	0.268322	0.00228627	0.0339351	0.0132941	-0.0188242	0.0245811
29	p57	0.0028259	-0.00255108	-0.0153236	0.0156019	0.0551938	-0.0361024
30	p53	-0.0783678	0.0257259	0.0203498	0.0180272	0.0430869	-0.0320897
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮

It appears to be Subj “p10”.

4.10.3 Item-related conditional means of random effects

Code

cm_Item = last(ranefinfo(m_prm2))
caterpillar!(Figure(; resolution=(800, 1200)), cm_Item; orderby=1)

Figure 4: Prediction intervals on item random effects for model m_prm2

5 Appendix

Code

versioninfo()

Julia Version 1.8.0
Commit 5544a0fab76 (2022-08-17 13:38 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin21.4.0)
  CPU: 12 × Intel(R) Core(TM) i9-8950HK CPU @ 2.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, skylake)
  Threads: 1 on 6 virtual cores